Method for dynamic state estimation of natural gas network considering dynamic characteristics of natural gas pipelines

ABSTRACT

Provided is a method for a dynamic state estimation of a natural gas network considering dynamic characteristics of natural gas pipelines. The method can obtain a result of the dynamic state estimation of the natural gas network by establishing an objective function of the dynamic state estimation of the natural gas network, a state quantity constraint of a compressor, a state quantity constraint of the natural gas pipeline and a topological constraint of the natural gas network, and using a Lagrange method or an interior point method to solve a state estimation model of the natural gas network. The method takes the topological constraint of the natural gas network into consideration, and employs a pipeline pressure constraint in a frequency domain to implement linearization of the pipeline pressure constraint, thereby obtain a real-time, reliable, consistent and complete dynamic operating state of the natural gas network.

TECHNICAL FIELD

The present disclosure relates to a method for a dynamic state estimation of a natural gas network considering dynamic characteristics of natural gas pipelines, belonging to the technical field of operation and control of an integrated energy system.

BACKGROUND

An integrated energy system has great advantages in terms of many aspects such as improving energy utilization efficiency, promoting new energy consumption, and reducing energy costs. It is a development trend of future energy systems. An Integrated Energy Management System (IEMS), which can regulate energy flows using information flows, is an intelligent decision-making “brain” that ensures a safe, economical, green, and highly efficient integrated energy system. The technology for estimating a state, as a basic module of the IEMS, is responsible for providing real-time, reliable, consistent, and complete operating state information, and thus it can provide the subsequent security analysis and optimization control with reliable operating data.

At present, researches on state estimations of a natural gas network are still in their infancy, not to mention state estimation technology that considers dynamic natural gas. Only some of the published literatures have proposed the dynamic state estimation methods based on Kalman filtering for a single natural gas pipeline. However, these methods fail to consider constraints of the natural gas network, and they require an initial state inside the pipeline to be known (generally assumed to be a steady state). In addition, a step of bad data identification is hard to be added to a format of an iterative solution of the Kalman filtering, which greatly limits its application. Thus, it is urgent to propose a method for a dynamic state estimation of a complex natural gas pipeline network, so as to provide sufficient data support for the operation and control of the integrated energy system.

SUMMARY

The present disclosure provides a method for a dynamic state estimation of a natural gas network considering dynamic characteristics of natural gas pipelines to obtain a real-time, reliable, consistent and complete operating state of the natural gas network, thereby overcoming the defects in the known state estimations of the natural gas network.

The method for the dynamic state estimation of the natural gas network considering the dynamic characteristic of the natural gas pipeline provided by the present disclosure includes:

step 1 of establishing a time-domain window and a frequency-domain window for the dynamic state estimation of the natural gas network, the step 1 including:

sub-step 1-1 of defining a time-domain window width as I^(t), where I^(t) is a positive integer, and a value of I^(t) is determined by a dispatcher of the natural gas network; defining a u-th sampling time point in the time-domain window as τ_(u)=τ−uΔt, u=0, 1, . . . , I^(t)−1, where τ represents a current time point of the natural gas network, and Δt represents a sampling interval of the natural gas network; defining a current time-domain window width as I^(t,e), where I^(t,e) is a positive integer, and a value of I^(t,e) is determined by the dispatcher of the natural gas network; and defining a historical time-domain window width as I^(t,h), where I^(t,h) is a positive integer, and a value of I^(t,h) is determined by the dispatcher of the natural gas network, wherein I^(t), I^(t,e) and I^(t,h) satisfy the following relational expression:

I ^(t) =I ^(t,e) +I ^(t,h); and

sub-step 1-2 of defining a frequency-domain window width as I^(f), where a value of I^(f) is determined by the dispatcher of the natural gas network; and defining a d-th frequency component in the frequency-domain window as ω_(d), d=0, 1, . . . , I^(f)−1, where ω_(d) is calculated by the following formula:

$\omega_{d} = \frac{d}{{I^{t} \cdot \Delta}\; t}$

step 2 of constructing a measurement vector for the dynamic state estimation of the natural gas network, the step 2 including:

sub-step 2-1 of acquiring, from a data acquisition and monitoring control system of the natural gas network, all operation data of the natural gas network at a sampling time point τ_(u) in the time-domain window where the current time point τ of the natural gas network belongs, wherein the all operation data of the natural gas network comprises: a measurement value z_(G) ₊ _(,u) ^(i) ^(p) of a natural gas flow at a head end of each pipeline in the natural gas network, and a measurement value z_(G) ⁻ _(,u) ^(i) ^(p) of a natural gas flow at a tail end of each pipeline in the natural gas network, where i_(p) represents a serial number of a pipeline in the natural gas network; a measurement value z_(G) ₊ _(,u) ^(i) ^(c) of a natural gas flow at a head end of each compressor, and a measurement value z_(G) ⁻ _(,u) ^(i) ^(c) of a natural gas flow at a tail end of each compressor, where i_(c) represents a serial number of a compressor; a pressure measurement value z_(pr,u) ^(i) ^(n) of each node of the natural gas network, where i_(n) represents a serial number of a node of the natural gas network; a measurement value z_(gs,u) ^(i) ^(s) of a natural gas flow of each natural gas source, where i_(s) represents a serial number of a natural gas source; and a measurement value z_(gl,u) ^(i) ^(l) of a natural gas flow of each natural gas load, where i_(l) represents a serial number of a natural gas load; and

sub-step 2-2 of constructing a measurement vector z_(u) for the dynamic state estimation of the natural gas network at the sampling time point τ_(u):

$Z_{u} = \begin{bmatrix} Z_{G^{+},u} \\ Z_{G^{-},u} \\ Z_{{pr},u} \\ Z_{{gs},u} \\ Z_{gl,u} \end{bmatrix}$

where z_(G) ₊ _(,u) represents a column vector consisting of all the measurement values z_(G) ₊ _(,u) ^(i) ^(p) of natural gas flows at head ends of respective pipelines in the natural gas network and all the measurement values z_(G) ₊ _(,u) ^(i) ^(c) of natural gas flows at head ends of respective compressors at the sampling time point τ_(u); z_(G) ⁻ _(,u) represents a column vector consisting of all the measurement values z_(G) ⁻ _(,u) ^(i) ^(p) of natural gas flows at tail ends of respective pipelines in the natural gas network and all the measurement values z_(G) ⁻ _(,u) ^(i) ^(c) of natural gas flows at tail ends of respective compressors at the sampling time point: τ_(u); z_(pr,u) represents a column vector consisting of all the pressure measurement values z_(pr,u) ^(i) ^(n) of respective nodes of the natural gas network at the sampling time point τ_(u); z_(gs,u) represents a column vector consisting of all the measurement values z_(gs,u) ^(i) ^(s) of natural gas flows of respective natural gas sources in the natural gas network at the sampling time point τ_(u); and z_(gl,u) represents a column vector consisting of all the measurement values z_(gl,u) ^(i) ^(l) of natural gas flows of respective natural gas loads in the natural gas network at the sampling time point τ_(u);

step 3 of constructing a state vector x_(u) for the dynamic state estimation of the natural gas network at the sampling time point τ_(u):

$x_{u} = \begin{bmatrix} x_{G^{+},u} \\ x_{G^{-},u} \\ x_{{pr},u} \\ x_{{gs},u} \\ x_{gl,u} \end{bmatrix}$

where x_(G) ⁻ _(,u) represents a column vector consisting of all natural gas flows G_(i) _(p) _(,u) ⁺ at the head ends of the respective pipelines in the natural gas network and all natural gas flows G_(i) _(c) _(,u) ⁺ at the head ends of the respective compressors at the sampling time point τ_(u); x_(G) ⁻ _(,u) represents a column vector consisting of all natural gas flows G_(i) _(p) _(,u) ⁺ at the tail ends of the respective pipelines in the natural gas network and all natural gas flows G_(i) _(c) _(,u) ⁻ at the tail ends of the respective compressors at the sampling time point τ_(u); x_(pr,u) represents a column vector consisting of all pressures h_(i) _(n) of the respective nodes of the natural gas network at the sampling time point τ_(u); x_(gs,u) represents a column vector consisting of all natural gas flows G_(i) _(s) _(,u) ^(gs) of the respective natural gas sources in the natural gas network at the sampling time point τ_(u); and x_(gl,u) represents a column vector consisting of all natural gas flows G_(i) _(l) _(,u) ^(gl) of the respective natural gas loads in the natural gas network at the sampling time point τ_(u);

step 4 of establishing, based on the measurement vector constructed in the step 2 and the state vector constructed in the step 3, an objective function of the dynamic state estimation of the natural gas network as follows:

min J=Σ _(u=0) ^(I) ^(t,e) ⁻¹{[z _(u) −x _(u)]W ⁻¹[z _(u) −x _(u)]^(T)}+Σ_(u=I) _(t,e) ^(I) ^(t) ⁻¹{[z _(u) −x _(u)]W ⁻¹δ^(u−I) ^(t,e) [z _(u) −x _(u)]^(T)}

where J represents an expression of the objective function; W represents a covariance matrix of a measurement error and is determined by the dispatcher of the natural gas network; a superscript T represents a matrix transpose; and δ represents a decay factor of a historical time window and is determined by the dispatcher of the natural gas network;

step 5 of establishing constraint conditions for the dynamic state estimation of the natural gas network, the step 5 including:

sub-step 5-1 of establishing constraints related to a flow and a pressure of the compressor in the natural gas network, the sub-step 5-1 including:

establishing a flow constraint of the head end and the tail end of a compressor:

G _(i) _(c) _(,u) ⁺ =G _(i) _(c) _(,u) ⁻ ,∀i _(c)ϵΩ_(c) ,∀u=0,1, . . . ,I ^(t)−1

where Ω_(c) represents a set of serial numbers of all the compressors in the natural gas network; and

establishing a pressure constraint at the head end and the tail end of the compressor, wherein, for the compressor with a constant tail end pressure, the pressure constraint of the head end and the tail end of the compressor is as follows:

h _(i) _(c) _(,u) ⁻ =h _(i) _(c) _(,con) ⁻ ,∀i _(c)ϵΩ_(c,1) ,∀u=0,1, . . . ,I ^(t)−1

where h_(i) _(c) _(,u) ⁻ represents a tail end pressure of a compressor i_(c) at the sampling time point τ_(u); h_(i) _(c) _(,con) ⁻ represents a set value of a tail end pressure of the compressor i_(c) and is a constant determined by the dispatcher of the natural gas network; and Ω_(c,1) represents a set of serial numbers of all the compressors with the constant tail end pressure in the natural gas network;

wherein, for a compressor with a constant compression ratio, the pressure constraint of the head end and the tail end of the compressor is as follows:

h _(i) _(c) _(,u) ⁻ =r _(i) _(c) _(,con) ·h _(i) _(c) _(,u) ⁻ ,∀i _(c)ϵΩ_(c,2) ,∀u=0,1, . . . ,I ^(t)−1

where h_(i) _(c) _(,u) ⁺ represents a head end pressure of the compressor i_(c) at the sampling time point τ_(u); r_(i) _(c) _(,con) represents a set value of a compression ratio of the compressor i_(c) and is a constant determined by the dispatcher of the natural gas network; and Ω_(c,2) represents a set of serial numbers of all the compressors with the constant compression ratio in the natural gas network; and

wherein, for a compressor with a constant pressure difference, the pressure constraint of the head end and the tail end of the compressor id as follows:

h _(i) _(c) _(,u) ⁻ −h _(i) _(c) _(,u) ⁺ =Δh _(i) _(c) _(,con) ,∀i _(c)ϵΩ_(c,3) ,∀u=0,1, . . . ,I ^(t)−1

where Δh_(i) _(c) _(,con) represents a set value of a pressure difference between the tail end and the head end of the compressor i_(c) and is a constant determined by the dispatcher of the natural gas network; and Φ_(c,3) represents a set of serial numbers of all the compressors with the constant pressure difference in the natural gas network;

sub-step 5-2 of establishing a flow constraint and a pressure constraint of the natural gas in the pipeline in the natural gas network, the sub-step 5-2 including:

establishing a two-port constraint of the pipeline in the natural gas network of each frequency component ω_(d) in the frequency-domain window:

${{\begin{bmatrix} h_{i_{p},d}^{-} \\ G_{i_{p},d}^{-} \end{bmatrix} = {\begin{bmatrix} A_{i_{p},d} & B_{i_{p},d} \\ C_{i_{p},d} & D_{i_{p},d} \end{bmatrix}\begin{bmatrix} h_{i_{p},d}^{+} \\ G_{i_{p},d}^{+} \end{bmatrix}}},{\forall{i_{p} \in \Omega_{p}}},{{\forall d} = 0},}{1,\ldots\mspace{14mu},{I^{f} - 1}}$

where h_(i) _(p) _(,d) ⁻ represents a value of a d-th component of a tail end pressure of a pipeline i_(p) in the natural gas network in a frequency-domain window I^(f), and h_(i) _(p) _(,d) ⁻ is a complex variable to be solved; h_(i) _(p) _(,d) ⁺ represents a value of a d-th frequency component of a head end pressure of the pipeline i_(p) in the natural gas network, and h_(i) _(p) _(,d) ⁺ is a complex variable to be solved; G_(i) _(p) _(,d) ⁻ represents a value of a d-th component of a natural gas flow at the tail end of the pipeline i_(p) in the natural gas network in the frequency-domain window I^(f), and G_(i) _(p) _(,d) ⁻ is a complex variable to be solved; G_(i) _(p) _(,d) ⁺ represents a value of a d-th component of a natural gas flow at the head end of the pipeline i_(p) in the natural gas network in the frequency-domain window, and G_(i) _(p) _(,d) ⁺ is a complex variable to be solved; and A_(i) _(p) _(,d), B_(i) _(p) _(,d), C_(i) _(p) _(,d) and D_(i) _(p) _(,d) represent two-port parameters of a d-th component of the pipeline i_(p) in the natural gas network in the frequency-domain window, and values of A_(i) _(p) _(,d), B_(i) _(p) _(,d), C_(i) _(p) _(,d) and D_(i) _(p) _(,d) are respectively expressed as:

$A_{i_{p},d} = {\left\lbrack {{\cosh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} - {b_{i_{p},d} \cdot {\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)}}} \right\rbrack \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}$ $B_{i_{p},d} = {{- \frac{z_{i_{p},d}}{\sqrt{\left( k_{i_{p}} \right)^{2} + {{4 \cdot Z_{i_{p},d}}Y_{i_{p},d}}}}}{{\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}}$ $C_{i_{p},d} = {{- \frac{Y_{i_{p},d}}{\sqrt{\left( k_{i_{p}} \right)^{2} + {{4 \cdot Z_{i_{p},d}}Y_{i_{p},d}}}}}{{\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}}$ $D_{i_{p},d} = {\left\lbrack {{\cosh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} + {b_{i_{p},d} \cdot {\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)}}} \right\rbrack \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}$

where l_(i) _(p) represents a length of the pipeline i_(p) in the natural gas network, k_(i) _(p) , a_(i) _(p) _(,d), b_(i) _(p) _(,d), Z_(i) _(p) _(,d) and Y_(i) _(p) _(,d) represent values of the d-th frequency component of pipeline parameters of the natural gas network, and values of k_(i) _(p) , a_(i) _(p) _(,d), b_(i) _(p) _(,d), Z_(i) _(p) _(,d) and Y_(i) _(p) _(,d) are respectively expressed as:

$k_{i_{p}} = {- \frac{{2gD_{i_{p}}\sin\;\theta_{i_{p}}} - {\lambda_{i_{p}}\left( v_{{base},i_{p}} \right)}^{2}}{2{RTD}_{i_{p}}}}$ $a_{i_{p},d} = {\frac{1}{2}\sqrt{\left( k_{i_{p},d} \right)^{2} + {4Z_{i_{p},d}Y_{i_{p},d}}}}$ $b_{i_{p},d} = \frac{k_{i_{p},d}}{\sqrt{\left( k_{i_{p},d} \right)^{2} + {4Z_{i_{p},d}Y_{i_{p},d}}}}$ Z_(i_(p), d) = R_(i_(p)) + jω_(d)L_(i_(p)) Y_(i_(p), d) = j ω_(d)C_(i_(p))

where g represents an acceleration of gravity; D_(i) _(p) represents an inner diameter of the pipeline i_(p) in the natural gas network; θ_(i) _(p) represents an angle of inclination of the pipeline i_(p) in the natural gas network; λ_(i) _(p) represents a friction coefficient of the pipeline i_(p) in the natural gas network; v_(base,i) _(p) represents a basic value of a flow velocity in the pipeline i_(p) in the natural gas network; R represents a gas constant of natural gas; T represents a temperature of the natural gas; j represents a complex number unit; and R_(i) _(p) , L_(i) _(p) and C_(i) _(p) are parameters of the natural gas network, and values of R_(i) _(p) , L_(i) _(p) and C_(i) _(p) are respectively expressed as:

R _(i) _(p) =λ_(i) _(p) v _(base,i) _(p) /(A _(i) _(p) D _(i) _(p) )

L _(i) _(p) =1/A _(i) _(p)

C _(i) _(p) =A _(i) _(p) /(RT)

where A_(i) _(p) represents a cross-sectional area of the pipeline i_(p) in the natural gas network;

establishing a time domain-frequency-domain mapping constraint of the natural gas flow at the head end of the pipeline of the natural gas network:

G _(i) _(p) _(,u) ⁺=Σ_(d=0) ^(I) ^(f) ⁻¹[Re(G _(i) _(p) _(,d) ⁺)·cos(θ_(d)−ω_(d) ·uΔt)−Im(G _(i) _(p) _(,d) ⁺)·sin(θ_(d)−ω_(d) ·uΔt)]

where Re( ) represents valuing a real part of a complex number; Im( ) represents valuing an imaginary part of the complex number; and θ_(d) represents a parameter calculated with ω_(d) as follows:

θ_(d) =I ^(f)·ω_(d)−ω_(d)

establishing a time domain-frequency-domain mapping constraint of the natural gas flow at the tail end of the pipeline of the natural gas network:

G _(i) _(p) _(,u) ⁻=Σ_(d=0) ^(I) ^(f) ⁻¹[Re(G _(i) _(p) _(,d) ⁻)·cos(θ_(d)−ω_(d) ·uΔt)−Im(G _(i) _(p) _(,d) ⁻)·sin(θ_(d)−ω_(d) ·uΔt)], and

establishing a time domain-frequency-domain mapping constraint of the node of the natural gas network:

h _(i) _(n) _(,u)=Σ_(d=0) ^(I) ^(f) ⁻¹[Re(h _(i) _(n) _(,d))·cos(θ_(d)−ω_(d) ·uΔt)−Im(h _(i) _(n) _(,d))·sin(θ_(d)−ω_(d) ·uΔt)]

where h_(i) _(n) _(,u) represents a value of a d-th component of a pressure of a node i_(n) in the frequency-domain window; and h_(i) _(n) _(,u) represents a complex variable to be solved; and

sub-step 5-3 of establishing a topological constraint of the natural gas network, the sub-step 5-3 including:

establishing a flow balance constraint of a node of the natural gas network:

∑_(i_(p) ∈ Ω_(p)^(+, i_(n)))G_(i_(p), u)⁺ − ∑_(i_(p) ∈ Ω_(p)^(−, i_(n)))G_(i_(p), u)⁻ + ∑_(i_(c) ∈ Ω_(c)^(+, i_(n)))G_(i_(c), u)⁺ − ∑_(i_(c) ∈ Ω_(c)^(−, i_(n)))G_(i_(c), u)⁻ + ∑_(i_(l) ∈ Ω_(l)^(−, i_(n)))G_(i_(l), u)^(gl) − ∑_(i_(s) ∈ Ω_(s)^(+, i_(n)))G_(i_(s), u)^(gs) = 0

where Ω_(p) ^(+,i) ^(n) represents a set of serial numbers of pipelines connected to the node i_(n) at head ends; Ω_(p) ^(−,i) ^(n) represents a set of serial numbers of pipelines connected to the node i_(n) at tail ends; Ω_(c) ^(+,i) ^(n) represents a set of serial numbers of compressors connected to the node i_(n) at head ends; Ω_(c) ^(−,i) ^(n) represents a set of serial numbers of compressors connected to the node i_(n) at tail ends; Ω_(s) ^(+,i) ^(n) represents a set of serial numbers of natural gas sources connected to the node i_(n); and Ω_(l) ^(−,i) ^(n) represents a set of serial numbers of natural gas loads connected to the node i_(n);

establishing constraints of a pipeline-compressor-node time-domain pressure relationship in the natural gas network:

h _(i) _(p) _(,u) ⁺ =h _(i) _(n) _(,u) ,∀i _(p)ϵΩ_(p) ^(+,i) ^(n)

h _(i) _(p) _(,u) ⁻ =h _(i) _(n) _(,u) ,∀i _(p)ϵΩ_(p) ^(−,i) ^(n)

h _(i) _(c) _(,u) ⁺ =h _(i) _(n) _(,u) ,∀i _(c)ϵΩ_(c) ^(−,i) ^(n)

h _(i) _(c) _(,u) ⁻ =h _(i) _(n) _(,u) ,∀i _(c)ϵΩ_(c) ^(−,i) ^(n) , and

establishing constraints of a pipeline-node frequency-domain pressure relationship in the natural gas network:

h _(i) _(p) _(,d) ⁺ =h _(i) _(n) _(,d) ,∀i _(p)ϵΩ_(p) ^(−,i) ^(n)

h _(i) _(p) _(,d) ⁻ =h _(i) _(n) _(,d) ,∀i _(p)ϵΩ_(p) ^(−,i) ^(n) , and

step 6 of forming a dynamic state estimation model of the natural gas network by using the objective function of the dynamic state estimation of the natural gas network established in the step 4 and the constraint conditions for the dynamic state estimation of the natural gas network established in the step 5; solving, by using a Lagrange method or an interior point method, the dynamic state estimation model of the natural gas network, to obtain the state vector x_(u), for the dynamic state estimation of the natural gas network at the sampling time point τ_(u); and performing the dynamic state estimation of the natural gas network by considering the dynamic characteristics of the natural gas pipelines.

The method provided by the present disclosure has the following advantages.

According to the present disclosure, the method for the dynamic state estimation of the natural gas network considering the dynamic characteristic of the natural gas pipeline can obtain a result of the dynamic state estimation of the natural gas network by establishing the objective function of the dynamic state estimation of the natural gas network, the state quantity constraint of the compressor, the state quantity constraint of the natural gas pipeline and the topological constraint of the natural gas network are established, and by using the Lagrange method or the interior point method to solve a state estimation model of the natural gas network. The method according to the present disclosure takes the topological constraint of the natural gas network into consideration, and employs a pipeline pressure constraint in a frequency domain to implement linearization of the pipeline pressure constraint, thereby obtaining a real-time, reliable, consistent and complete dynamic operating state of the natural gas network.

DESCRIPTION OF EMBODIMENTS

A method for a dynamic state estimation of a natural gas network considering dynamic characteristics of natural gas pipelines provided by the present disclosure includes the following steps (1) to (5).

(1) A time-domain window and a frequency-domain window for the dynamic state estimation of the natural gas network are established. The step (1) includes the following steps (1-1) to (1-2).

(1-1) A time-domain window width is defined as I^(t), where I^(t) is a positive integer, and a value of I^(t) is determined by a dispatcher of the natural gas network. A u-th sampling time point in the time-domain window is defined as τ_(u)=τ−uΔt, u=0, 1, . . . , I^(t)−1, where τ represents a current time point of the natural gas network, and Δt represents a sampling interval of the natural gas network. A current time-domain window width is defined as I^(t,e), where I^(t,e) is a positive integer, and a value of I^(t,e) is determined by the dispatcher of the natural gas network. A historical time-domain window width is defined as I^(t,h), where I^(t,h) is a positive integer. A value of I^(t,h) is determined by the dispatcher of the natural gas network. I^(t), I^(t,e) and I^(t,h) satisfy the following relational expression:

I ^(t) =I ^(t,e) +I ^(t,h).

(1-2) A frequency-domain window width is defined as I^(f), where a value of I^(f) is determined by the dispatcher of the natural gas network. A d-th frequency component in the frequency-domain window is defined as ω_(d), d=0, 1, . . . , I^(f)−1, where ω_(d) is calculated by the following formula:

${\omega_{d} = \frac{d}{{I^{t} \cdot \Delta}\; t}}.$

(2) A measurement vector for the dynamic state estimation of the natural gas network is constructed. The step (2) includes the following steps (2-1) to (2-1).

(2-1) All operation data of the natural gas network at a sampling time point τ_(u) in the time-domain window where the current time point T of the natural gas network belongs is acquired from a data acquisition and monitoring and control system of the natural gas network. The all operation data of the natural gas network includes: a measurement value z_(G) ₊ _(,u) ^(i) ^(p) of a natural gas flow at a head end of each pipeline in the natural gas network, and a measurement value z_(G) ⁻ _(,u) ^(i) ^(p) of a natural gas flow at a tail end of each pipeline in the natural gas network, where i_(p) represents a serial number of a pipeline in the natural gas network; a measurement value z_(G) ₊ _(,u) ^(i) ^(c) of a natural gas flow at a head end of each compressor, and a measurement value z_(G) ⁻ _(,u) ^(i) ^(c) of a natural gas flow at a tail end of each compressor, where i_(c) represents a serial number of a compressor; a pressure measurement value z_(pr,u) ^(i) ^(n) of each node of the natural gas network, where i_(n) represents a serial number of a node of the natural gas network; a measurement value z_(gs,u) ^(i) ^(s) of a natural gas flow of each natural gas source, where i_(s) represents a serial number of a natural gas source; and a measurement value z_(gl,u) ^(i) ^(l) of a natural gas flow of each natural gas load, where i_(l) represents a serial number of a natural gas load.

(2-2) A measurement vector z_(u) for the dynamic state estimation of the natural gas network at each sampling time point τ_(u) is constructed:

${z_{u} = \begin{bmatrix} Z_{G^{+},u} \\ Z_{G^{-},u} \\ Z_{{pr},u} \\ Z_{{gs},u} \\ Z_{{gl},u} \end{bmatrix}},$

where z_(G) ₊ _(,u) represents a column vector consisting of all the measurement values z_(G) ₊ _(,u) ^(i) ^(p) of natural gas flows at head ends of respective pipelines in the natural gas network and all the measurement values z_(G) ₊ _(,u) ^(i) ^(c) of natural gas flows at head ends of respective compressors at the sampling time point τ_(u); z_(G) ⁻ _(,u) represents a column vector consisting of all the measurement values z_(G) ⁻ _(,u) ^(i) ^(p) of natural gas flows at tail ends of respective pipelines in the natural gas network and all the measurement values z_(G) ⁻ _(,u) ^(i) ^(c) of natural gas flows at tail ends of respective compressors at the sampling time point τ_(u); z_(pr,u) represents a column vector consisting of all the pressure measurement values z_(pr,u) ^(i) ^(n) of respective nodes of the natural gas network at the sampling time point τ_(u); z_(gs,u) represents a column vector consisting of all the measurement values z_(gs,u) ^(i) ^(s) of natural gas flows of respective natural gas sources in the natural gas network at the sampling time point τ_(u); and z_(gl,u) represents a column vector consisting of all the measurement values z_(gl,u) ^(i) ^(l) of natural gas flows of respective natural gas loads in the natural gas network at the sampling time point τ_(u).

(3) A state vector x_(u) for the dynamic state estimation of the natural gas network at each sampling time point τ_(u) is constructed:

${x_{u} = \begin{bmatrix} x_{G^{+},u} \\ x_{G^{-},u} \\ x_{{pr},u} \\ x_{{gs},u} \\ x_{{gl},u} \end{bmatrix}},$

where x_(G) ⁻ _(,u) represents a column vector consisting of all natural gas flows G_(i) _(p) _(,u) ⁺ at the head ends of the respective pipelines in the natural gas network and all natural gas flows G_(i) _(c) _(,u) ⁺ at the head ends of the respective compressors at the sampling time point τ_(u); x_(G) ⁻ _(,u) represents a column vector consisting of all natural gas flows G_(i) _(p) _(,u) ⁻ at the tail ends of the respective pipelines in the natural gas network and all natural gas flows G_(i) _(c) _(,u) ⁻ at the tail ends of the respective compressors at the sampling time point τ_(u), x_(pr,u) represents a column vector consisting of all pressures h_(i) _(n) of the respective nodes of the natural gas network at the sampling time point τ_(u); x_(gs,u) represents a column vector consisting of all natural gas flows G_(i) _(s) _(,u) ^(gs) of the respective natural gas sources in the natural gas network at the sampling time point τ_(u); and x_(gl,u) represents a column vector consisting of all natural gas flows G_(i) _(l) _(,u) ^(gl) of the respective natural gas loads in the natural gas network at the sampling time point τ_(u).

(4) An objective function of the dynamic state estimation of the natural gas network is established based on the measurement vector constructed in step (2) and the state vector constructed in step (3):

min J=Σ _(u=0) ^(I) ^(t,e) ⁻¹{[z _(u) −x _(u)]W ⁻¹[z _(u) −x _(u)]^(T)}+Σ_(u=I) _(t,e) ^(I) ^(t) ⁻¹{[z _(u) −x _(u)]W ⁻¹δ^(u−I) ^(t,e) [z _(u) −x _(u)]^(T)}

where J represents an expression of the objective function; W represents a covariance matrix of a measurement error and is determined by the dispatcher of the natural gas network; a superscript T represents a matrix transpose; and δ represents a decay factor of a historical time window and is determined by the dispatcher of the natural gas network.

(5) Constraint conditions for the dynamic state estimation of the natural gas network are established. The step (5) includes steps (5-1) to (5-3).

(5-1) Constraints related to a flow and a pressure of a compressor in the natural gas network are established. The step (5-1) includes steps (5-1-1) to (5-1-2).

(5-1-1) A flow constraint at a head end and a tail end of a compressor is established:

G _(i) _(c) _(,u) ⁺ =G _(i) _(c) _(,u) ⁻ ,∀i _(c)ϵΩ_(c) ,∀u=0,1, . . . ,I ^(t)−1

where Ω_(c) represents a set of serial numbers of respective compressors in the natural gas network.

(5-1-2) A pressure constraint at a head end and a tail end of a compressor is established.

For a compressor with a constant tail end pressure, the pressure constraint at the head end and the tail end of the compressor being as follows:

h _(i) _(c) _(,u) ⁻ =h _(i) _(c) _(,con) ⁻ ,∀i _(c)ϵΩ_(c,1) ,∀u=0,1, . . . ,I ^(t)−1

where h_(i) _(c) _(,u) ⁻ represents a tail end pressure of a compressor i_(c) at the sampling time point τ_(u), h_(i) _(c) _(,con) ⁻ represents a set value of the tail end pressure of the compressor i_(c), is a constant, and is determined by the dispatcher of the natural gas network, and Ω_(c,1) represents a set of serial numbers of respective compressors with the constant tail end pressure in the natural gas network.

For a compressor with a constant compression ratio, the pressure constraint of the head end and the tail end of the compressor is as follows:

h _(i) _(c) _(,u) ⁻ =r _(i) _(c) _(,con) ·h _(i) _(c) _(,u) ⁻ ,∀i _(c)ϵΩ_(c,2) ,∀u=0,1, . . . ,I ^(t)−1

where h_(i) _(c) _(,u) ⁺ represents a head end pressure of the compressor i_(c) at the sampling time point τ_(u); r_(i) _(c) _(,con) represents a set value of a compression ratio of the compressor i_(c) and is a constant determined by the dispatcher of the natural gas network; and Ω_(c,2) represents a set of serial numbers of all the compressors with the constant compression ratio in the natural gas network.

For a compressor with a constant pressure difference, the pressure constraint of the head end and the tail end of the compressor is as follows:

h _(i) _(c) _(,u) ⁻ −h _(i) _(c) _(,u) ⁺ =Δh _(i) _(c) _(,con) ,∀i _(c)ϵΩ_(c,3) ,∀u=0,1, . . . ,I ^(t)−1

where Δh_(i) _(c) _(,con) represents a set value of a pressure difference between the tail end and the head end of the compressor i_(c) and is a constant determined by the dispatcher of the natural gas network; and Φ_(c,3) represents a set of serial numbers of all the compressors with the constant pressure difference in the natural gas network.

(5-2) A flow constraint and a pressure constraint of natural gas in a pipeline in the natural gas network are established. The step (5-2) includes steps (5-2-1) to (5-2-5).

(5-2-1) A two-port constraint of the pipeline in the natural gas network of each frequency component ω_(d) in the frequency-domain window is established:

${\begin{bmatrix} h_{i_{p},d}^{-} \\ G_{i_{p},d}^{-} \end{bmatrix} = {\begin{bmatrix} A_{i_{p},d} & B_{i_{p},d} \\ C_{i_{p},d} & D_{i_{p},d} \end{bmatrix}\begin{bmatrix} h_{i_{p},d}^{+} \\ G_{i_{p},d}^{+} \end{bmatrix}}},{{\text{∀}i_{p}} \in \Omega_{p}},{{\text{∀}d} = 0},1,\ldots\mspace{14mu},{I^{f} - 1},$

where h_(i) _(p) _(,d) ⁻ represents a value of a d-th component of a tail end pressure of a pipeline i_(p) in the natural gas network in a frequency-domain window I^(f), and h_(i) _(p) _(,d) ⁻ is a complex variable to be solved; h_(i) _(p) _(,d) ⁺ represents a value of a d-th frequency component of a head end pressure of the pipeline i_(p) in the natural gas network, and h_(i) _(p) _(,d) ⁺ is a complex variable to be solved; G_(i) _(p) _(,d) ⁻ represents a value of a d-th component of a natural gas flow at the tail end of the pipeline i_(p) in the natural gas network in the frequency-domain window I^(f), and G_(i) _(p) _(,d) ⁻ is a complex variable to be solved; G_(i) _(p) _(,d) ⁺ represents a value of a d-th component of a natural gas flow at the head end of the pipeline i_(p) in the natural gas network in the frequency-domain window, and G_(i) _(p) _(,d) ⁻ is a complex variable to be solved; and A_(i) _(p) _(,d), B_(i) _(p) _(,d), C_(i) _(p) _(,d) and D_(i) _(p) _(,d) represent two-port parameters of a d-th component of the pipeline i_(p) in the natural gas network in the frequency-domain window, and values of A_(i) _(p) _(,d), B_(i) _(p) _(,d), C_(i) _(p) _(,d) and D_(i) _(p) _(,d) are respectively expressed as:

$A_{i_{p},d} = {\left\lbrack {{\cosh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} - {b_{i_{p},d} \cdot {\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)}}} \right\rbrack \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}$ $B_{i_{p},d} = {{- \frac{Z_{i_{p},d}}{\sqrt{\left( k_{i_{p}} \right)^{2} + {{4 \cdot Z_{i_{p},d}}Y_{i_{p},d}}}}}{{\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}}$ $C_{i_{p},d} = {{- \frac{Y_{i_{p},d}}{\sqrt{\left( k_{i_{p}} \right)^{2} + {{4 \cdot Z_{i_{p},d}}Y_{i_{p},d}}}}}{{\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}}$ ${D_{i_{p},d} = {\left\lbrack {{\cosh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} - {b_{i_{p},d} \cdot {\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)}}} \right\rbrack \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}},$

where l_(i) _(p) represents a length of the pipeline i_(p) in the natural gas network, k_(i) _(p) , a_(i) _(p) _(,d), b_(i) _(p) _(,d), Z_(i) _(p) _(,d) and Y_(i) _(p) _(,d) represent values of the d-th frequency component of pipeline parameters of the natural gas network, and values of k_(i) _(p) , a_(i) _(p) _(,d), b_(i) _(p) _(,d), Z_(i) _(p) _(,d) and Y_(i) _(p) _(,d) are respectively expressed as:

$k_{i_{p}} = {- \frac{{2gD_{i_{p}}\sin\;\theta_{i_{p}}} - {\lambda_{i_{p}}\left( v_{{base},i_{p}} \right)}^{2}}{2{RTD}_{i_{p}}}}$ $a_{i_{p},d} = {\frac{1}{2}\sqrt{\left( k_{i_{p},d} \right)^{2} + {4Z_{i_{p},d}Y_{i_{p},d}}}}$ $b_{i_{p},d} = \frac{k_{i_{p},d}}{\sqrt{\left( k_{i_{p},d} \right)^{2} + {4Z_{i_{p},d}Y_{i_{p},d}}}}$ Z_(i_(p), d) = R_(i_(p)) + jω_(d)L_(i_(p)) Y_(i_(p), d) = j ω_(d)C_(i_(p))

where g represents an acceleration of gravity; D_(i) _(p) represents an inner diameter of the pipeline i_(p) in the natural gas network; θ_(i) _(p) represents an angle of inclination of the pipeline i_(p) in the natural gas network; λ_(i) _(p) represents a friction coefficient of the pipeline i_(p) in the natural gas network; v_(base,i) _(p) represents a basic value of a flow velocity in the pipeline i_(p) in the natural gas network; R represents a gas constant of natural gas; T represents a temperature of the natural gas; j represents a complex number unit; and R_(i) _(p) , L_(i) _(p) and C_(i) _(p) are parameters of the natural gas network, and values of R_(i) _(p) , L_(i) _(p) and C_(i) _(p) are respectively expressed as:

R _(i) _(p) =λ_(i) _(p) v _(base,i) _(p) /(A _(i) _(p) D _(i) _(p) )

L _(i) _(p) =1/A _(i) _(p)

C _(i) _(p) =A _(i) _(p) /(RT)

where A_(i) _(p) represents a cross-sectional area of the pipeline i_(p) in the natural gas network.

(5-2-3) A time domain-frequency-domain mapping constraint of the natural gas flow at the head end of the pipeline of the natural gas network is established:

G _(i) _(p) _(,u) ⁺=Σ_(d=0) ^(I) ^(f) ⁻¹[Re(G _(i) _(p) _(,d) ⁺)·cos(θ_(d)−ω_(d) ·uΔt)−Im(G _(i) _(p) _(,d) ⁺)·sin(θ_(d)−ω_(d) ·uΔt)],

where Re( ) represents valuing a real part of a complex number; Im( ) represents valuing an imaginary part of the complex number; and θ_(d) represents a parameter calculated with ω_(d) as follows:

θ_(d) =I ^(f)·ω_(d)−ω_(d).

(5-2-4) A time domain-frequency-domain mapping constraint of the natural gas flow at the tail end of the pipeline of the natural gas network is established:

G _(i) _(p) _(,u) ⁻=Σ_(d=0) ^(I) ^(f) ⁻¹[Re(G _(i) _(p) _(,d) ⁻)·cos(θ_(d)−ω_(d) ·uΔt)−Im(G _(i) _(p) _(,d) ⁻)·sin(θ_(d)−ω_(d) ·uΔt)].

(5-2-5) A time domain-frequency-domain mapping constraint of nodes of the natural gas network is established:

h _(i) _(n) _(,u)=Σ_(d=0) ^(I) ^(f) ⁻¹[Re(h _(i) _(n) _(,d))·cos(θ_(d)−ω_(d) ·uΔt)−Im(h _(i) _(b) _(,d))·sin(θ_(d)−ω_(d) ·uΔt)],

where h_(i) _(n) _(,u) represents a value of a d-th component of a pressure of a node i_(n) in the frequency-domain window; and h_(i) _(n) _(,u) represents a complex variable to be solved.

(5-3) A topological constraint of the natural gas network is established. The step (5-3) includes steps (5-3-1) to (5-3-3).

(5-3-1) A flow balance constraint of a node of the natural gas network is established:

∑_(i_(p) ∈ Ω_(p)^(+, i_(n)))G_(i_(p), u)⁺ − ∑_(i_(p) ∈ Ω_(p)^(−, i_(n)))G_(i_(p), u)⁻ + ∑_(i_(c) ∈ Ω_(c)^(+, i_(n)))G_(i_(c), u)⁺ − ∑_(i_(c) ∈ Ω_(c)^(−, i_(n)))G_(i_(c), u)⁻ + ∑_(i_(l) ∈ Ω_(l)^(−, i_(n)))G_(i_(l), u)^(gl) − ∑_(i_(s) ∈ Ω_(s)^(+, i_(n)))G_(i_(s), u)^(gs) = 0,

where Ω_(p) ^(+,i) ^(n) represents a set of serial numbers of pipelines connected to the node i_(n) at head ends; Ω_(p) ^(−,i) ^(n) represents a set of serial numbers of pipelines connected to the node i_(n) at tail ends; Ω_(c) ^(+,i) ^(n) represents a set of serial numbers of compressors connected to the node i_(n) at head ends; Ω_(c) ^(−,i) ^(n) represents a set of serial numbers of compressors connected to the node i_(n) at tail ends; Ω_(s) ^(+,i) ^(n) represents a set of serial numbers of natural gas sources connected to the node i_(n); and Ω_(l) ^(−,i) ^(n) represents a set of serial numbers of natural gas loads connected to the node i_(n).

(5-3-2) Constraints of a pipeline-compressor-node time-domain pressure relationship in the natural gas network are established:

h _(i) _(p) _(,u) ⁺ =h _(i) _(n) _(,u) ,∀i _(p)ϵΩ_(p) ^(+,i) ^(n)

h _(i) _(p) _(,u) ⁻ =h _(i) _(n) _(,u) ,∀i _(p)ϵΩ_(p) ^(−,i) ^(n)

h _(i) _(c) _(,u) ⁺ =h _(i) _(n) _(,u) ,∀i _(c)ϵΩ_(c) ^(−,i) ^(n)

h _(i) _(c) _(,u) ⁻ =h _(i) _(n) _(,u) ,∀i _(c)ϵΩ_(c) ^(−,i) ^(n) .

(5-3-3) Constraints of a pipeline-node frequency-domain pressure relationship in the natural gas network are established:

h _(i) _(p) _(,d) ⁺ =h _(i) _(n) _(,d) ,∀i _(p)ϵΩ_(p) ^(−,i) ^(n)

h _(i) _(p) _(,d) ⁻ =h _(i) _(n) _(,d) ,∀i _(p)ϵΩ_(p) ^(−,i) ^(n) .

(6) A dynamic state estimation model of the natural gas network is formed by using the objective function of the dynamic state estimation of the natural gas network established in step (4) and the constraint conditions for the dynamic state estimation of the natural gas network established in step (5). The dynamic state estimation model of the natural gas network is solved by using a Lagrange method or an interior point method, to obtain the state vector x_(u) for the dynamic state estimation of the natural gas network at each sampling time point τ_(u). Consequently, the dynamic state estimation of the natural gas network considering the dynamic characteristic of the natural gas pipeline is implemented.

In the embodiments of the present disclosure, commercial software Gurobi or Cplex is used to solve the dynamic state estimation model of the natural gas network. 

What is claimed is:
 1. A method for a dynamic state estimation of a natural gas network considering dynamic characteristics of natural gas pipelines, the method comprising: step 1 of establishing a time-domain window and a frequency-domain window for the dynamic state estimation of the natural gas network, the step 1 comprising: sub-step 1-1 of defining a time-domain window width as I^(t), where I^(t) is a positive integer, and a value of I^(t) is determined by a dispatcher of the natural gas network; defining a u-th sampling time point in the time-domain window as τ_(u)=τ−uΔt, u=0, 1, . . . , I^(t)−1, where τ represents a current time point of the natural gas network, and Δt represents a sampling interval of the natural gas network; defining a current time-domain window width as I^(t,e), where I^(t,e) is a positive integer, and a value of I^(t,e) is determined by the dispatcher of the natural gas network; and defining a historical time-domain window width as I^(t,h), where I^(t,h) is a positive integer, and a value of I^(t,h) is determined by the dispatcher of the natural gas network, wherein I^(t), I^(t,e) and I^(t,h) satisfy the following relational expression: I ^(t) =I ^(t,e) +I ^(t,h); and sub-step 1-2 of defining a frequency-domain window width as I^(f), where a value of I^(f) is determined by the dispatcher of the natural gas network; and defining a d-th frequency component in the frequency-domain window as ω_(d), d=0, 1, . . . , I^(f)−1, where ω_(d) is calculated by the following formula: ${\omega_{d} = \frac{d}{{I^{t} \cdot \Delta}\; t}}.$ step 2 of constructing a measurement vector for the dynamic state estimation of the natural gas network, the step 2 comprising: sub-step 2-1 of acquiring, from a data acquisition and monitoring control system of the natural gas network, all operation data of the natural gas network at a sampling time point τ_(u) in the time-domain window where the current time point r of the natural gas network belongs, wherein the all operation data of the natural gas network comprises: a measurement value z_(G) ₊ _(,u) ^(i) ^(p) of a natural gas flow at a head end of each pipeline in the natural gas network, and a measurement value z_(G) ⁻ _(,u) ^(i) ^(p) of a natural gas flow at a tail end of each pipeline in the natural gas network, where i_(p) represents a serial number of a pipeline in the natural gas network; a measurement value z_(G) ₊ _(,u) ^(i) ^(c) of a natural gas flow at a head end of each compressor, and a measurement value z_(G) ⁻ _(,u) ^(i) ^(c) of a natural gas flow at a tail end of each compressor, where i_(c) in represents a serial number of a compressor; a pressure measurement value z_(pr,u) ^(i) ^(n) of each node of the natural gas network, where i_(n) represents a serial number of a node of the natural gas network; a measurement value z_(gs,u) ^(i) ^(s) of a natural gas flow of each natural gas source, where i_(s) represents a serial number of a natural gas source; and a measurement value z_(gl,u) ^(i) ^(l) of a natural gas flow of each natural gas load, where i_(l) represents a serial number of a natural gas load; and sub-step 2-2 of constructing a measurement vector z_(u) for the dynamic state estimation of the natural gas network at the sampling time point τ_(u): ${z_{u} = \begin{bmatrix} Z_{G^{+},u} \\ Z_{G^{-},u} \\ Z_{{pr},u} \\ Z_{{gs},u} \\ Z_{{gl},u} \end{bmatrix}},$ where z_(G) ₊ _(,u) represents a column vector consisting of all the measurement values z_(G) ₊ _(,u) ^(i) ^(p) of natural gas flows at head ends of respective pipelines in the natural gas network and all the measurement values z_(G) ₊ _(,u) ^(i) ^(c) of natural gas flows at head ends of respective compressors at the sampling time point τ_(u); z_(G) ⁻ _(,u) represents a column vector consisting of all the measurement values z_(G) ⁻ _(,u) ^(i) ^(p) of natural gas flows at tail ends of respective pipelines in the natural gas network and all the measurement values z_(G) ⁻ _(,u) ^(i) ^(c) of natural gas flows at tail ends of respective compressors at the sampling time point τ_(u); z_(pr,u) represents a column vector consisting of all the pressure measurement values z_(pr,u) ^(i) ^(n) of respective nodes of the natural gas network at the sampling time point τ_(u); z_(gs,u) represents a column vector consisting of all the measurement values z_(gs,u) ^(i) ^(s) of natural gas flows of respective natural gas sources in the natural gas network at the sampling time point τ_(u); and z_(gl,u) represents a column vector consisting of all the measurement values z_(gl,u) ^(i) ^(l) of natural gas flows of respective natural gas loads in the natural gas network at the sampling time point τ_(u); step 3 of constructing a state vector x_(u) for the dynamic state estimation of the natural gas network at the sampling time point τ_(u): ${x_{u} = \begin{bmatrix} x_{G^{+},u} \\ x_{G^{-},u} \\ x_{{pr},u} \\ x_{{gs},u} \\ x_{{gl},u} \end{bmatrix}},$ where x_(G) ₊ _(,u) represents a column vector consisting of all natural gas flows G_(i) _(p) _(,u) ⁺ at the head ends of the respective pipelines in the natural gas network and all natural gas flows G_(i) _(c) _(,u) ⁺ at the head ends of the respective compressors at the sampling time point τ_(u); x_(G) ⁻ _(,u) represents a column vector consisting of all natural gas flows G_(i) _(p) _(,u) ⁻ at the tail ends of the respective pipelines in the natural gas network and all natural gas flows G_(i) _(c) _(,u) ⁻ at the tail ends of the respective compressors at the sampling time point τ_(u); x_(pr,u) represents a column vector consisting of all pressures h_(i) _(n) of the respective nodes of the natural gas network at the sampling time point τ_(u); x_(gs,u) represents a column vector consisting of all natural gas flows G_(i) _(s) _(,u) ^(gs) of the respective natural gas sources in the natural gas network at the sampling time point τ_(u); and x_(gl,u) represents a column vector consisting of all natural gas flows G_(i) _(l) _(,u) ^(gl) of the respective natural gas loads in the natural gas network at the sampling time point τ_(u); step 4 of establishing, based on the measurement vector constructed in the step 2 and the state vector constructed in the step 3, an objective function of the dynamic state estimation of the natural gas network as follows: min J=Σ _(u=0) ^(I) ^(t,e) ⁻¹{[z _(u) −x _(u)]W ⁻¹[z _(u) −x _(u)]^(T)}+Σ_(u=I) _(t,e) ^(I) ^(t) ⁻¹{[z _(u) −x _(u)]W ⁻¹δ^(u−I) ^(t,e) [z _(u) −x _(u)]^(T)}, where J represents an expression of the objective function; W represents a covariance matrix of a measurement error and is determined by the dispatcher of the natural gas network; a superscript T represents a matrix transpose; and δ represents a decay factor of a historical time window and is determined by the dispatcher of the natural gas network; step 5 of establishing constraint conditions for the dynamic state estimation of the natural gas network, the step 5 comprising: sub-step 5-1 of establishing constraints related to a flow and a pressure of the compressor in the natural gas network, the sub-step 5-1 comprising: establishing a flow constraint of the head end and the tail end of a compressor: G _(i) _(c) _(,u) ⁺ =G _(i) _(c) _(,u) ⁻ ,∀i _(c)ϵΩ_(c) ,∀u=0,1, . . . ,I ^(t)−1 where Ω_(c) represents a set of serial numbers of all the compressors in the natural gas network; and establishing a pressure constraint at the head end and the tail end of the compressor, wherein, for the compressor with a constant tail end pressure, the pressure constraint of the head end and the tail end of the compressor is as follows: h _(i) _(c) _(,u) ⁻ =h _(i) _(c) _(,con) ⁻ ,∀i _(c)ϵΩ_(c,1) ,∀u=0,1, . . . ,I ^(t)−1 where h_(i) _(c) _(,u) ⁻ represents a tail end pressure of a compressor i_(c) at the sampling time point τ_(u); h_(i) _(c) _(,con) ⁻ represents a set value of a tail end pressure of the compressor i_(c) and is a constant determined by the dispatcher of the natural gas network; and Ω_(c,1) represents a set of serial numbers of all the compressors with the constant tail end pressure in the natural gas network; wherein, for a compressor with a constant compression ratio, the pressure constraint of the head end and the tail end of the compressor is as follows: h _(i) _(c) _(,u) ⁻ =r _(i) _(c) _(,con) ·h _(i) _(c) _(,u) ⁺ ,∀i _(c)ϵΩ_(c,2) ,∀u=0,1, . . . ,I ^(t)−1 where h_(i) _(c) _(,u) ⁻ represents a head end pressure of the compressor i_(c) at the sampling time point τ_(u); r_(i) _(c) _(,con) represents a set value of a compression ratio of the compressor i_(c) and is a constant determined by the dispatcher of the natural gas network; and Ω_(c,2) represents a set of serial numbers of all the compressors with the constant compression ratio in the natural gas network; and wherein, for a compressor with a constant pressure difference, the pressure constraint of the head end and the tail end of the compressor id as follows: h _(i) _(c) _(,u) ⁻ −h _(i) _(c) _(,u) ⁺ =Δh _(i) _(c) _(,con) ,∀i _(c)ϵΩ_(c,3) ,∀u=0,1, . . . ,I ^(t)−1 where Δh_(i) _(c) _(,con) represents a set value of a pressure difference between the tail end and the head end of the compressor i_(c) and is a constant determined by the dispatcher of the natural gas network; and Φ_(c,3) represents a set of serial numbers of all the compressors with the constant pressure difference in the natural gas network; sub-step 5-2 of establishing a flow constraint and a pressure constraint of the natural gas in the pipeline in the natural gas network, the sub-step 5-2 comprising: establishing a two-port constraint of the pipeline in the natural gas network of each frequency component ω_(d) in the frequency-domain window: ${\begin{bmatrix} h_{i_{p},d}^{-} \\ G_{i_{p},d}^{-} \end{bmatrix} = {\begin{bmatrix} A_{i_{p},d} & B_{i_{p},d} \\ C_{i_{p},d} & D_{i_{p},d} \end{bmatrix}\begin{bmatrix} h_{i_{p},d}^{+} \\ G_{i_{p},d}^{+} \end{bmatrix}}},{{\text{∀}i_{p}} \in \Omega_{p}},{{\text{∀}d} = 0},1,\ldots\mspace{14mu},{I^{f} - 1},$ where h_(i) _(p) _(,d) ⁻ represents a value of a d-th component of a tail end pressure of a pipeline i_(p) in the natural gas network in a frequency-domain window I^(f), and h_(i) _(p) _(,d) ⁻ is a complex variable to be solved; h_(i) _(p) _(,d) ⁺ represents a value of a d-th frequency component of a head end pressure of the pipeline i_(p) in the natural gas network, and h_(i) _(p) _(,d) ⁺ is a complex variable to be solved; G_(i) _(p) _(,d) ⁻ represents a value of a d-th component of a natural gas flow at the tail end of the pipeline i_(p) in the natural gas network in the frequency-domain window I^(f), and G_(i) _(p) _(,d) ⁻ is a complex variable to be solved; G_(i) _(p) _(,d) ⁺ represents a value of a d-th component of a natural gas flow at the head end of the pipeline i_(p) in the natural gas network in the frequency-domain window, and G_(i) _(p) _(,d) ⁺ is a complex variable to be solved; and A_(i) _(p) _(,d), B_(i) _(p) _(,d), C_(i) _(p) _(,d) and D_(i) _(p) _(,d) represent two-port parameters of a d-th component of the pipeline i_(p) in the natural gas network in the frequency-domain window, and values of A_(i) _(p) _(,d), B_(i) _(p) _(,d), C_(i) _(p) _(,d) and D_(i) _(p) _(,d) are respectively expressed as: $A_{i_{p},d} = {\left\lbrack {{\cosh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} - {b_{i_{p},d} \cdot {\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)}}} \right\rbrack \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}$ $B_{i_{p},d} = {{- \frac{Z_{i_{p},d}}{\sqrt{\left( k_{i_{p}} \right)^{2} + {{4 \cdot Z_{i_{p},d}}Y_{i_{p},d}}}}}{{\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}}$ $C_{i_{p},d} = {{- \frac{Y_{i_{p},d}}{\sqrt{\left( k_{i_{p}} \right)^{2} + {{4 \cdot Z_{i_{p},d}}Y_{i_{p},d}}}}}{{\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}}$ ${D_{i_{p},d} = {\left\lbrack {{\cosh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)} - {b_{i_{p},d} \cdot {\sinh\left( {l_{i_{p}} \cdot a_{i_{p},d}} \right)}}} \right\rbrack \cdot e^{- \frac{k_{i_{p}} \cdot l_{i_{p}}}{2}}}},$ where l_(i) _(p) represents a length of the pipeline i_(p) in the natural gas network, k_(i) _(p) , a_(i) _(p) _(,d), b_(i) _(p) _(,d), Z_(i) _(p) _(,d) and Y_(i) _(p) _(,d) represent values of the d-th frequency component of pipeline parameters of the natural gas network, and values of k_(i) _(p) , a_(i) _(p) _(,d), b_(i) _(p) _(,d), Z_(i) _(p) _(,d) and Y_(i) _(p) _(,d) are respectively expressed as: $k_{i_{p}} = {- \frac{{2gD_{i_{p}}\sin\;\theta_{i_{p}}} - {\lambda_{i_{p}}\left( v_{{base},i_{p}} \right)}^{2}}{2{RTD}_{i_{p}}}}$ $a_{i_{p},d} = {\frac{1}{2}\sqrt{\left( k_{i_{p},d} \right)^{2} + {4Z_{i_{p},d}Y_{i_{p},d}}}}$ $b_{i_{p},d} = \frac{k_{i_{p},d}}{\sqrt{\left( k_{i_{p},d} \right)^{2} + {4Z_{i_{p},d}Y_{i_{p},d}}}}$ Z_(i_(p), d) = R_(i_(p)) + jω_(d)L_(i_(p)) Y_(i_(p), d) = j ω_(d)C_(i_(p)) where g represents an acceleration of gravity; D_(i) _(p) represents an inner diameter of the pipeline i_(p) in the natural gas network; θ_(i) _(p) represents an angle of inclination of the pipeline i_(p) in the natural gas network; λ_(i) _(p) represents a friction coefficient of the pipeline i_(p) in the natural gas network; v_(base,i) _(p) represents a basic value of a flow velocity in the pipeline i_(p) in the natural gas network; R represents a gas constant of natural gas; T represents a temperature of the natural gas; j represents a complex number unit; and R_(i) _(p) , L_(i) _(p) and C_(i) _(p) are parameters of the natural gas network, and values of R_(i) _(p) , L_(i) _(p) and C_(i) _(p) are respectively expressed as: R _(i) _(p) =λ_(i) _(p) v _(base,i) _(p) /(A _(i) _(p) D _(i) _(p) ) L _(i) _(p) =1/A _(i) _(p) C _(i) _(p) =A _(i) _(p) /(RT) where A_(i) _(p) represents a cross-sectional area of the pipeline i_(p) in the natural gas network; establishing a time domain-frequency-domain mapping constraint of the natural gas flow at the head end of the pipeline of the natural gas network: G _(i) _(p) _(,u) ⁺=Σ_(d=0) ^(I) ^(f) ⁻¹[Re(G _(i) _(p) _(,d) ⁺)·cos(θ_(d)−ω_(d) ·uΔt)−Im(G _(i) _(p) _(,d) ⁺)·sin(θ_(d)−ω_(d) ·uΔt)], where Re( ) represents valuing a real part of a complex number; Im( ) represents valuing an imaginary part of the complex number; and θ_(d) represents a parameter calculated with ω_(d) as follows: θ_(d) =I ^(f)·ω_(d)−ω_(d); establishing a time domain-frequency-domain mapping constraint of the natural gas flow at the tail end of the pipeline of the natural gas network: G _(i) _(p) _(,u) ⁻=Σ_(d=0) ^(I) ^(f) ⁻¹[Re(G _(i) _(p) _(,d) ⁻)·cos(θ_(d)−ω_(d) ·uΔt)−Im(G _(i) _(p) _(,d) ⁻)·sin(θ_(d)−ω_(d) ·uΔt)]; and establishing a time domain-frequency-domain mapping constraint of the node of the natural gas network: $h_{i_{n},u} = {\sum\limits_{d = 0}^{I^{f} - 1}\left\lbrack {{{{Re}\left( h_{i_{n},d} \right)} \cdot {\cos\left( {\theta_{d} - {{\omega_{d} \cdot u}\;\Delta\; t}} \right)}} - {{{Im}\left( h_{i_{n},d} \right)} \cdot {\sin\left( {\theta_{d} - {{\omega_{d} \cdot u}\;\Delta\; t}} \right\rbrack}}} \right.}$ where h_(i) _(n) _(,u) represents a value of a d-th component of a pressure of a node i_(n) in the frequency-domain window; and h_(i) _(n) _(,u) represents a complex variable to be solved; and sub-step 5-3 of establishing a topological constraint of the natural gas network, the sub-step 5-3 comprising: establishing a flow balance constraint of a node of the natural gas network: ${{\sum\limits_{i_{p} \in \Omega_{p}^{+ {,i_{n}}}}G_{i_{p},u}^{+}} - {\sum\limits_{i_{p} \in \Omega_{p}^{- {,i_{n}}}}G_{i_{p},u}^{-}} + {\sum\limits_{i_{c} \in \Omega_{c}^{+ {,i_{n}}}}G_{i_{c},u}^{+}} - {\sum\limits_{i_{c} \in \Omega_{p}^{- {,i_{n}}}}G_{i_{p},u}^{-}} + {\sum\limits_{i_{l} \in \Omega_{l}^{- {,i_{n}}}}G_{i_{l},u}^{gl}} - {\sum\limits_{i_{s} \in \Omega_{s}^{+ {,i_{n}}}}G_{i_{s},u}^{gs}}} = 0$ where Ω_(p) ^(+,i) ^(n) represents a set of serial numbers of pipelines connected to the node i_(n) at head ends; Ω_(p) ^(−,i) ^(n) represents a set of serial numbers of pipelines connected to the node i_(n) at tail ends; Ω_(c) ^(+,i) ^(n) represents a set of serial numbers of compressors connected to the node i_(n) at head ends; Ω_(c) ^(−,i) ^(n) represents a set of serial numbers of compressors connected to the node i_(n) at tail ends; Ω_(s) ^(+,i) ^(n) represents a set of serial numbers of natural gas sources connected to the node i_(n); and Ω_(l) ^(+,i) ^(n) represents a set of serial numbers of natural gas loads connected to the node i_(n); establishing constraints of a pipeline-compressor-node time-domain pressure relationship in the natural gas network: h _(i) _(p) _(,u) ⁺ =h _(i) _(n) _(,u) ,∀i _(p)ϵΩ_(p) ^(+,i) ^(n) h _(i) _(p) _(,u) ⁻ =h _(i) _(n) _(,u) ,∀i _(p)ϵΩ_(p) ^(−,i) ^(n) h _(i) _(c) _(,u) ⁺ =h _(i) _(n) _(,u) ,∀i _(c)ϵΩ_(c) ^(−,i) ^(n) h _(i) _(c) _(,u) ⁻ =h _(i) _(n) _(,u) ,∀i _(c)ϵΩ_(c) ^(−,i) ^(n) ; and establishing constraints of a pipeline-node frequency-domain pressure relationship in the natural gas network: h _(i) _(p) _(,d) ⁺ =h _(i) _(n) _(,d) ,∀i _(p)ϵΩ_(p) ^(−,i) ^(n) h _(i) _(p) _(,d) ⁻ =h _(i) _(n) _(,d) ,∀i _(p)ϵΩ_(p) ^(−,i) ^(n) ; and step 6 of forming a dynamic state estimation model of the natural gas network by using the objective function of the dynamic state estimation of the natural gas network established in the step 4 and the constraint conditions for the dynamic state estimation of the natural gas network established in the step 5; solving, by using a Lagrange method or an interior point method, the dynamic state estimation model of the natural gas network, to obtain the state vector x_(u) for the dynamic state estimation of the natural gas network at the sampling time point τ_(u); and performing the dynamic state estimation of the natural gas network by considering the dynamic characteristics of the natural gas pipelines. 